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1 . Review 

a. Physical problems and test computations 

Our research has concentrated on the numerical ‘Solution" of a single, 
mixed, nonlinear equation with prescribed boundary data. The governing 
equation has the general form 

{C(y) + 4> x > <t> xx - = 0 (1) 

For the caustic problem, which we examined first, C(y) has the simple form 

C(y) - y. (2) 

and the boundary data are prescribed in accordance with the asymptotic 
analog of (1), viz., 

^xx ~ ^yy 

These boundary data include an incoming signal, 

<P X = " P y " 1/4 f (p)> alon 9 q = - 1 , (3b) 

i 

where F(p) gives the shape of the incoming signal, p characterizes its strength, 
and p,q = x + yy^^, respectively. 

For the two-dimensional transonic flow problems, which we are now investi- 
gating, C(y) is the constant transonic parameter, - K. The boundary data here 
are the far-field behavior of the solution and the tangency condition on the 
airfoil . 

We have developed a second-order numerical procedure for "solving 11 (1) 
and a shock-fitting scheme to treat the discontinuities that appear in the 
solution. Numerical results for this (the caustic) problem are codified in 
the appended paper, which will be delivered at the First International Confe- 
rence on Computational Method in Nonlinear Mechanics, September 23-25, 1974. 
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In the following sections, we outline the research progress we have made 
to date and further work now underway. These include a brief review of the 
second order scheme, the shock fitting technique, comparison of the results 
of various schemes, and a preliminary test computation of a transonic flow 
problem. 


b. Second-order scheme 

The second-order scheme which is used in the present computation is 
similar to the first-order scheme of Murman and Cole (1971), i.e., in elliptic 
regions (C(y) + <b < 0) a central difference approximation is chosen to 

A 

represent the x- ‘and the y-derivatives : 

= <Vu- -Vi,p /2Ax ’ 

♦*x “' ( *1+l.j ■ 2 *1.j + *1-1.3 ,/to2 ’ (4a) 

*yy = ( *i . j+l ‘ 2< h ,j + * 1 , j-l )/Ay2 ' 


In hyperbolic regions (C(y) + 4> > 0) the x-derivatives are replaced by one- 

A 

sided differences to correctly represent the domain of dependence. 


= (2<f, i.j " ~ 24> i-2,j + ■ 

*xx = (2 *|,j - 5 Vl,j + 4 *i-2,j 1 < *'i-3,j )/AxZ • 

After substituting (4a) and (4b) into (1), we obtain a set of nonlinear 
difference equations of the form 

A 1 <] + l + A 2 *1J + A 3 <]-l + A 5 = V* 5 = 0 • 


(4b) 


(5) 


The coefficients An, are shown in Table 1 of the appended paper. Equation (5) 
is flux conservative ; we solve it by the Newton's method: 

(<f> n+1 - $ n ) aFj/3^ = - FjU) . 


( 6 ) 



In our last report we pointed out that the matrix 8F./9<f>. loses diagonal 
dominance in hyperbolic compressive regions, which causes poor convergence or 
instability in the iterative procedure when the initial guess is not accurate. 
We remedied this by introducing an artificial viscosity, v<j> , such that (1) 

aaA 

becomes 

{C(y) + 4> X H XX “ <Pyy ~ W> xxx * (?) 

The diagonal term of (6) becomes 

3F a /8<,, 0 = ^2 + ^2 {C(y) + *X } + to ^xx + V/Ax3 ; (8) 


v is zero when 


&{y) + <f> ) + 
Ax^ ■ 


Ax 



> 0 . 


(9a) 


When this is not the case, we choose 

v = kAx2 {k [C(y) + ^x 3 + v) * (9b) 

where k > 1 

With this modification, equation (6) has diagonal dominance everywhere in the 

computational field. Test computations show that the present scheme provides a 

converged solution for a wide range of signal strengths (y = 0.05 to 0.25) 

with the linear solution to (3a) as an initial guess. After the numerical value 

of <p. . becomes more accurate (usually about. 20 computations*), we can set 
i > J 

v - 0 and continue the computational procedure without difficulty until the 
"solution" converges. The results for our second-order calculation for signal 
strength of y = 0.05 are displayed in Figure 3 of the appendix. 


A computation is one complete calculation of the flow field. 



c. Shock fitting scheme for the caustic problem 

Usual finite difference schemes give poor representations of flow fields 
with shock waves embedded in them. This is the consequence of replacing deri- 
vatives by difference approximations in the whole computation region regard- 
less of the discontinuous behavior of the solution. A correct "solution" can 
only be obtained by treating shocks as a discontinuity. Moretti (1972,1973) 
has emphasized the importance of shock fitting for flows containing shock 
waves in several of his papers. Without shock fitting a large number of grid 
points and considerable computing time are required to achieve a given accuracy 
In the appendix we demonstrate the computational advantage of shock-fitting 
scheme by studying the nonlinear acoustic behavior of a shock wave near a 
caustic using different numerical schemes. The difficulties associated with 
shock fitting include determining the location of the shock and simplifying 
the numerical program so that it can be easily used. We have developed a 
successful scheme to solve the caustic problem; for other applications, our 
results require proper general ization. 

We assume a discontinuous signal with prescribed shape and strength given 
by (3b). The initial position of the signal (shock) is known (Figure 1 of 
the appendix) and the successive position of the shock along each computation 
column x = x. is determined by 

(dy/dx) s = t (C(y) + 1 (cf> x + J x )} ' 1/2 . (10) 

The flow properties ahead of the shock are obtained by the characteristic 
relations 

(c(y) + <f> x ) 1/2 d(f> x = + , (ID 

and the flow properties right behind the shock are determined by one of 
equations (11) and the jump condition 

($ x ■ <f> x ) 2 Wy) + \ (<J> X + 5> x )} = ($ v - . 


02 ) 
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Table 2 of the appendix details the procedure used to compute the position of 

shock and determine <f>, <J> , <p at the shock. 

x y 

Difference equations are then constructed by considering the shock point 
as a grid point. The entire difference equations including both shock, and 
shock free, regions are listed in Table 1 of the appendix. At every computa- 
tion point, the type dependent coefficient, C(y) + <j> is evaluated by a central 

A 

difference approximation. At the same time, the position of the shock is de- 
tected such that the proper difference equation can be chosen. The difference 
equation in shock region is not in conservative form, since such a form is not 
possible. Internal sources or sinks due to the nonconservative approximation 
are localized and can be considered negligible as we have avoided differencing 
across sharp gradients. 

In the region where the downstream condition of the shock becomes subsonic, 
equation (11) is no longer valid and we have to evaluate the flow property, 
e.g., $ , by a one-sided difference approximation, and calculate $ by (12). 

V 

As the incident shock reaches the sonic line, i.e., C(y) + <j> x = 0, we 
assume a reflected wave is formed. The initial strength of the reflected shock 
is obtained by using backward differences to approximate properties ahead of 
the shock and using forward differences to approximate properties behind the 
shock. With such an approximation, the reflected shock is quite weak. We have 
tried local symmetrical shocks at the triple point as proposed in the last 
progress report, but these do not seem to have the correct properties as we 
have not obtained a converged solution using them. A more general triple-point 
model is now under investigation. 

Computer drawn plots of 4> x /u at different y are shown in Figure 5 of the 


appendix. 



2. Comparisons of various schemes 


In Figure 2 of the appendix the linear solution of equation (3) for 
U = 0.05, and 6 = 20 (Gill and Seebass, 1974} is reproduced. This serves as 
an initial guess for our second-order scheme. Figure 4 of the appendix shows 
the first-order numerical results of Seebass, Murman, and Krupp (1971). First- 
order truncation errors diffuse the shock waves, especially near triple point 
where the incident shock and the reflected shock can not be clearly distin- 
guished. Figure 3 of the appendix is the result obtained by the present 
second-order scheme. The shock waves become sharper and the amplitude of the 
shock jump increases. However, in the region where dispersive errors dominate, 
i.e., |C(y) + <j> | > |<f> | , an unpleasant "wiggle” appears downstream of the 

shock. Such phenomenon can be eliminated by fitting a shock wave that satis- 
fies the correct jump relations. Figure 5 of the appendix shows the result 
of the shock-fitting scheme after 20 computations from the converged second- 
order "solution". The dispersive "wiggle" vanishes after the first computa- 
tion with the shock-fitting scheme. Table 3 of the appendix shows the dif- 
ferences in computing time and the rate of convergence for different schemes. 
Compared with the first-order and the second-order schemes, computing time 
per computation with the shock-fitting scheme is doubled, but the total 
computing time to bring the result to the same final accuracy is appreciably 
reduced. The fast rate of convergence and the correct shock jump are the 
main advantages of our shock fitting scheme. 

3. Transonic flow over two dimensional airfoils 

a. Governing equation 

The first-order approximation for transonic flow over a thin airfoil is 
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{K - V *xx + ♦nn = 0- O 3 ) 

where K is the transonic similarity parameter, and is defined as 

K = 0 - n£)/(2r ] t (14) 

and $ is the nondimensional perturbed velocity potential, which is related 
to the perturbed velocities u', v', by 

(21- M 2 ) 1/3 

V = - 2/3 ( U ‘/U,) (15a) 

T 

= 7 (V/Uj. (15b) 

Here r is the thickness ratio, = (y + l)/2, x, n are related to the 
physical coordinates by 

X = x/c, n = (2r 1 T m 2 ) 1/3 y/c, (16) 

where c is the chord length. 

The boundary conditions are then on y = t Y(x), or n -* 0 

<f> n = V(x), (17a) 

and at infinity, 

. *X ’ + 0 

for < 1 . (17b) 


For > 1 far-field boundary data are calculated by the method of character- 

istics. For example, at boundary point B in the following sketch, we have 

<*x - K s > 3/2 d *x = - d V 

or 

2/3 {(* - K ) 3/Z - (4 - K ) 3/2 } = - * R + 4 

x s g X s m ri M 

along the characteristic Z^. 
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Both <P Xoo and .6 are zero on C 2 , thus 

*nB “ f «- k ) 3/Z - (* x - U 3 b /2 >- (170 
Using (13), we can calculate by a one- 
sided difference approximation, i.e.. 



- 1.0 


1 1.0 X 



(17d) 


where <t> is computed by a backward difference approximation. 

X B 

For a finite computation region, the far field boundary data for < 1 
are obtained by using the formula derived by Murman arsd Cole (1971), 


1 f 1 

<j> « — . 2 f 
2tt s/K -1 


x-C 


(x-£) +K.n 


2 Y(£)dC 


x +K.n -» 


If dSdn 


i 


(18) 


The first term in (18) is the traditional Prandtl-Glauert solution for the 
linear asymptotic analog of (13); and the second term is the nonlinear cor- 
rection due to the perturbed velocity <j> . 

On the airfoil $ is approximated by a reflected boundary condition, 

i.e., 

% = i U 2 - 2 *1 + V = ^2 { 2 <*2 -* 1 >- ^ \ j , 09 ) 

where subscript "l" refers to the point on the body. 

The same numerical procedure discussed in section 1 is followed. A test 
computation for a 6% thick parabolic, nonlifting airfoil at various Mach 
numbers has been carried out. For subsonic flow over the airfoil, the far 
field solution is corrected after every 5 forward computations. At the same 



time the computation procedure is reversed, i.e., from downstream to upstream 
to bring the effect of the downstream boundary condition to the whole flow 
field as soon as practicable. Converged solutions have been obtained using 
our second-order scheme for M OT = 0.806, 0.861, 0.909. Figure 1 shows the 
present results, the results of first-order computation, and the results of 
Murman (1974). There is a major discrepancy in the shock wave position 
between our results and those calculated by Murman. This results in a large 
disagreement in the wave drag:. 

Murman and Cole (1974): Cp = 0.0315 

Present: 0.0125 

Knetchtel : 0.00835 

This discrepancy could be attributed to our use of a normal shock on the 
shock-fitting scheme, but initial calculations which allow for a floating 
shock of general shape, locate the shock at x/c = 0.875 t 0.025. This dis- 
agreement needs to be resolved and we are in the process of refining our cal- 
culations in order to do so. 

b. Embedded shock 

For supercritical flow, i.e., subsonic flow with embedded supersonic 
region, a shock wave is formed near the trailing edge through coalescence of 
waves originating on the body and reflected from the sonic line. Numerical 
results for the 6% parabolic arc airfoil show a rapid compression region near 
the trailing edge for M<» > 0.85. Correct representation of embedded shock 
requires a shock-fitting scheme. We assume the shock originates in far-field 
with infinitesimal strength, where the position of the shock is determined 
by the following criteria: 

(K - <j> x ) c = 0 + , and (K - <f> x ) b = O' ; (20) 

the subscript ”c" refers to the quantity evaluated by central difference 



• Mur man's shock point operator 
+ Experiment* NASA TN D-15 
© First order solutions 
— Present results 


vs. x for 6% arc nonlifting airfoil 
at M * 0.909, n = 0. 

oo 1 
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approximation, and "b" the backward difference approximation. This criteria 
is equivalent to insisting that the shock form in a supersonic, compression 
region. Once the starting point of the shock is determined, we can then trace 
out the shock position by the procedure of section lc. Our preliminary compu- 
tations have been limited to fitting a normal shock throughout the flow field. 
This is not a true representation of the shock that appears in the flow, but 
it provides some qualitative features of the flow field. Figure 2 shows the 
velocity variation at different n for = 0.909. An expansion occurs behind 
the shock; this probably corresponds to the Zierep (1958) singularity which 
does not show up in many numerical results. A floating shock fitting scheme 
is now under investigation. 


c. Wave drag 

Transonic wave drag can either be calculated by a contour integration of 
the pressure around the airfoil, or by the entropy production due to shock 
waves. Murman and Cole (1974) have pointed out. that due to the leading edge 
singularity, the embedded shock waves, and the fast variation of pressure near 
trailing edge, an accurate drag is hard to obtain by contour integration around 
the body. The alternative method is to compute the entropy production due to 
shock waves, and then to apply Oswatitsch's drag formula to determine the 
wave drag, 

oo 

D = 1 im p T / (s - s ) dy 

n OO OO J ' 00' J 

X-KO -00 


= P« T co/ , shock dy - 
shocks 


( 21 ) 


We know that (13) is the correct first-order approximation for transonic 
flow over a thin airfoil in the domain where shock waves and other discontinuities 
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are excluded. In the sketch, S represents the 
shock, W the wake, and B the airfoil. D is 
the domain where (13) applies. A weak solution 
of (13) provides a proper jump in pressure and 
velocity across S, which can be used to relate 
to the entropy change to the pressure jump: 

2 

As /c = 1 V (P/P - l) 3 (22) 

V 12 y 

The pressure jump (P/P) - 1 can be replaced by 
the velocity jump through 



or 


P/P - 1 = UP/P J/ (P/P J - U 

= {1 + y/2 M 2 (- 2 u'/Uj) {1 + y/2 M 2 (- 2 u'/Uj}- 1 
= - Y "l (u'/U ra - u'/Uj 


2 3 „4 
y M 


. a , 3 . 3 g , ^ 3 ^ T m /v ~ 

(P/P-l ) 3 = - Y 3 (u'/U„ - U'/UJ 3 = - 2 f - -- (4> x - *j 3 . 


After substituting (22) and (24) into (21), we have 


(23) 

(24) 



- p U 2 M 2 t 5 ' 3 
00 00 00 

12 (21^ M^) 1/3 


/ ($ x - ) 3 dn . 

shocks * * 


(25) 


Once the shock jump <j> x - <f> is accurately calculated by a shock-fitting scheme, 
the drag can easily be evaluated from (25). A different approach in deriving 
(25) can be found in Murman and Cole (1974). 
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4. Applications and further research 

Most two-dimensional problems with shock waves can be studied by the 
present shock-fitting scheme. For flow problems with multiple shocks, special 
care has to be paid to the region where shock waves merge or intersect each 
other. For three dimensional problems, the scheme becomes more complicated; 
we intend to try a simple, finite span wing problem, with shock waves 
determined by the converged second-order solution. We also intend to compare 
the result of our scheme with Moretti's three-dimensional shock-fitting result 
(1973). 

At present, we are improving our calculation of flow over a two-dimensional 
airfoil using the shock-fitting scheme we have developed. For the caustic 
problem we already studied, numerical computations with a much finer mesh near 
the triple point will be tested. The goal of the present research is to 
generalize the present numerical procedure to be useful for most two-dimensional 
and some simple three-dimensional flow problems. 
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Abstract 


This paper discusses the procedures we have developed to treat a can- 
onical problem involving a mixed nonlinear equation with boundary data 
that imply a discontinuous solution. This equation arises in various phys- 
ical contexts and is basic to the description of the nonlinear acoustic 
behavior of a shock wave near a caustic. The numerical scheme developed 
is of second order, treats discontinuities as such by applying the appro- 
priate jump conditions across them, and eliminates the numerical dissipa- 
tion and dispersion associated with large gradients. Our results are com- 
pared with the results of a first-order scheme and with those of a second- 
order scheme we have developed. The algorithm used here can easily be 
generalised to more complicated problems, including transonic flows with 
imbedded shocks. 


Introduction 

The computation of mixed, e.g., transonic, flows has been investigat- 
ed extensively in the past decade. Recent surveys of the numerical proce- 
dures used can be found in Nieuwland and Spee (1973), and Yoshihara (1972). 
The numerical treatment of such flows when shock waves are present has not 
been developed satisfactorily. Recently, Murman (1974) improved previous 
relaxation procedures by introducing a "shock point operator" to the dif- 
ference equations. His scheme notes the local character of the flow, and 
provides a relaxation scheme that insures that the calculation is fully 
conservative. However, due to first-order truncation errors, shock waves 
are smeared out, and consequently shock wave geometry can not be accurate- 
ly predicted. We present a numerical procedure for solving mixed equat- 
ions with second-order numerical accuracy by treating discontinuities as 
such. Moretti (1969 , 1972, 1973) has pursued a similar course in tack- 
ling related problems. 

Governing Equation and Boundary Conditions 

(y + $)<i> -(f) =0 

Y x' xx yy 


We consider 


0 ) 



where <J> may be thought of as a (perturbed) velocity potential, with bound- 
ary data prescribed in accordance with the nonlinear generalization of 
properly posed problems for 


yd> -<6 =0. 

J ^xx ^yy 


(2a) 


Numerical boundary data are determined from the solution to (2a) for an 
incoming signal with 


<l> x ^ - uy 1/4 F(p) 


for q 


(2b) 


where F(p) gives the shape of an incoming signal, u characterizes its 
strength, and p, q = x t 2/3 y 3 ' 2 , respectively. Equation (1) arises in 
various physical contexts; one of these is discussed in some detail in See- 
bass (1971). A discontinuous signal with strength \x = 0.05 was chosen for 
the present study with 

F(p) = H(p) - H(p + 6), (2c) 

where H(p) is the Heaviside unit function, and <5 was taken to be 20. This 
problem is sketched in figure 1 for. the domain considered here. 

The characteristic directions and the corresponding compatibility 
relations for (1) are _ 1/P 

dy/dx = t (y + d> x ) _ ' » (3) 

(y + *x> 1/2 * x = t d<j) y . (4) 

If discontinuities are present in the solution of (1), then they must have 
the directions * mo 

dy/dx = + {y + f (4> x + ^ x )> 7 (5). 

and across them 

(?x ” + \ ^x + = ^y " ^y^ ( 6 ) 

where, e.g.» - $ is the jump in <£ across the discontinuity. 

Solutions to the linear problem (2) may be calculated with any pre- 
cision desired (Gill and Seebass, 1974). Values of for fixed y are 
displayed in figure 2. The results provide an initial guess of the solu- 
tion to (1), as well as the boundary data. The computation is carried out 
in the region of figure 1. At points on the boundary where y + $ < 0, <jj 

is prescribed; at appropriate points on the boundary where y + <J> x > 0»4>and 
f are prescribed. On certain portions of the boundary, no data* are pre- 
scribed because the solution is determined uniquely without them. 


A first-order numerical "solution" to (1) was obtained by Seebass, Wur- 
man, and Krupp (1971) with an implicit, backward difference approximation 
to x-derivatives chosen for the grid points that lie in the hyperbolic re- 
gion. The schema is unconditionally stable and the numerical calculations 
converge. However, the "solution" does not give a satisfactory represents- 



tion of the discontinuities. We have developed a modified second-order 
scheme that solves (1), the smeared “discontinuities" obtained are con- 
siderably thinner than those obtained from the first-order scheme. One 
drawback of the second-order scheme is that in certain regions dispersive 
errors dominate and an unpleasant "wiggle" appears on the "downstream" 
side of the discontinui ty. • Using these second-order results for initial 
conditions, we then proceed with a second-order "shock-fitting" scheme 
that treats the discontinuities as such in order to satisfy the jump 
conditions to second-order. 


Numerical Procedure 


Second-Order Scheme 

The difference equations for (1) are of the form 

A 1 ♦i?Jll. + fl 2<J +A 3+i!!jll +A 5 = 0 


(7) 


or'Fj(<|>).= 0, where the index "i" refers to grid points in the x-direction, 
“j" J to the grid points in the y-direction, and the superscript "n+l" to 
the number of iterations of the entire region. The coefficients, An, are 
listed in table 1. Equation (7) is linear in elliptic domains, and non- 
linear in hyperbolic domains and can be solved by Newton's method, i.e.. 


4> n+1 - 4> n = -'OF.M/a*.} -1 F {<#.). 

J k J 


H - f'NC /x \ /<u i ' tr / jl .\ 

The difference approximation displayed in table 1 has the truncation 


error 


I A>(2 <>xx*xxx + W A * 2(y + ♦jPW 


1 9 

tW Ay (I) , for y + A > 0, 

12 J *yyyy y x 


- Tr Ax** 6 $ - Ax 2 (y + d> )d> +■ JL- Ay*~d> 

6 v xx\xxx 12 v x /y xxxx 12 " y 


(9) 


, for y + A < 0. 

yyyy J 


Special care' has to be taken when Newton's method is applied to hyp- 
erbolic domains. For hyperbolic equations the numerical error will not 
decay unless a proper scheme has been used. In the present problem the 
diagonal term of the tri-diagonal matrix 3F.{$)/3$. for hyperbolic domains 

. ~ < J K 


is of the form 


2 


Diag(j) = + 

A y c Ax^ Ax 


2 to* s n+1 + A n+1 > 

3 + ^i-2, j } 


~~ + _2_ ( + . ) + 1 . 

A 2 A 2 x' Ax y xx 

Ay A x^ 


(10a) 


In hyperbolic compressive regions the matrix loses diagonal dominance when 

7 (y + *x } + ^ *xx 
Ax* 

3 


< 0 . 


(10b) 



This can lead either to poor convergence on to instability of Newton's 
method. Thus an amendment is made by adding an artificial term k<J) , to 
(1), such that the difference equation remains diagonally dominant x/x in 
the iterative procedure. The value of k is of second order, and is deter- 
mined by 


Ax' 


+ *x) + It *xx 


> 0 


(VI) 


Ax' 


With this modification, a stable second-order numerical solution is obtain- 
ed. Computer- drawn plots of -<j> /p at constant y are shown in figure 3. Whai 
these results are compared witn the first-order computation shown in figure 
4, it can be seen that the second-order scheme provides sharper and thinner 
"shocks 21 . However, as mentioned, in regions where | y + 4> x l > l ( t ) xx l > "wig- 
gles" appear. . 


Shock-Fitting Technique 

Moretti (1972) has emphasized the importance of treating discontinui- 
ties as such for flows containing shock waves. He calls this procedure 
“shock fitting". Without shock fitting, a large number of grid points and 
considerable computing time are required to achieve a given accuracy. The 
advantage of shock fitting is clearly demonstrated by the present study. 
With shock fitting, we can determine both the shock position and the shock 
pressure rise with reasonable accuracy. Moreover, effects of numerical dis- 
sipation and dispersion are reduced to a minimum. 

We assume the computation procedure has reached station "i", i * e . , 
x = x-, as shown in the sketch of table 2; the upstream conditions are than 
all kJiown, and the properties of the downstream shock point "b" can be cal- 
culated by using the characteristic relations along bd, be, bf, and the 
jump condition (6). At point "c", where the shock intersects vertical grid 
line x - x., the value of <j> is calculated by direct integration of d$ alorg 
be. We thJn construct the difference approximation to x, and y derivative 
by using the shock points b and c instead of the regular grid points, e.g., 
h and k.. Again, an implicit scheme is used when the equation is hyperbolic 
and a central -difference scheme when it is elliptic. During the computa- 
tion the position of the shock is determined and the quantity y + <J> is 
computed at each grid point so that the proper difference equations' can 
be selected. 


As the incident discontinuity approaches the line where y + 6 - 0 it 
grows in strength until the flow right behind the incident wave Becomes 
sonic. At that point wa have assumed that a reflected wave is formed. The 
initial strength of the reflected wave can be obtained either by construct- 
ing local symmetrical discontinuities that meet at this point (a triple 
point), or by using backward differences to approximate properties ahead 
of the shock and using forward differences to approximate properties behird 
the shock. The second method is easier to use and the present study is 
limited to it. The strength of the discontinuity increases rapidly and ap- 
proaches its final value within 5 to 10. iterations. Computer-drawn plots 
of - (h v/ /u at constant y are shown in figure 5. i ha reflected discontinuity 
is weak, but this may be due, in part, to the procedure we invoked at the 



triple point. However, the reflected wave is not weak and is poorly repre- 
sented by the second-order results without shock fitting. 

Conclusion 

The numerical scheme outlined here offers a reliable method of comput- 
ing solutions to the mixed nonlinear equations with discontinuities. Compa- 
rison of the graphical results for different schemes shows the present 
method provides quantitatively superior results for an equal investment in 
computational time. 

Stability criteria for shock-fitting procedures derive from the same 
arguments used for the first-order and second-order implicit scheme. The 
rate of convergence may be studied by examining the maximum error of <t> for 
each successive calculation along column x. which requires the maximum 
number of iterations to compute, i.e.. 

Maximum error = !??? , |l - .j along x. . 

J —» ,u I , J l , J Mil 

Table 3 compares the computation time and the maximum error for the first- 
order, the second-order , and the shock-fitting schemes. It can be seen that 
for 60 iterations both the first-order and the second-order schemes have 
approximately the same rate of convergence. The initial guess was the lin- 
ear solution for both calculations. We can not use the linear solution 'as 
initial data with shock-fitting because the linear solution is too poor an 
initial guess. Our computation used the results from the second-order 
scheme as initial data. It took only twenty computations to reduce the max- 
imum error to 1%. The fast rate of convergence probably derives from the 
accuracy of the second-order solution away from the discontinuities, but it 
also indicates the efficiency of our procedure. The most undesirable feat- 
ure of the shock-fitting scheme is that the program becomes complicated . 
with bookkeeping. However, for two-dimensional problems, even those with 
multiple discontinuities, the present scheme seems easy to apply. For three- 
dimensional problems the difficulties are more substantial. 
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Computer drawn plots of - $ /y at 
constant y for the linear solution. 
Increments in y are 0.25. 


Computer drawn plots of - $ x /\i at 
constancy for the second-order ! 
calculation. Increments in y are 0.25. 
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Computer drawn plots of - <J> x /u at Computer drawn plots of - 4> x /p at 

constant y for a first-order calcu- constant y with shock fitting, 

lation. Increments in y are 0.25. Increments in y are 0.25. 
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Table K Formulation of Difference Equations 


A 1 + U+1 + A 2 <j + A 3 * lt]-l + Vs + A 5 * 0 


Case 1, Shock-free region: 


/ «r A. 1 /A n+ 1 OA n+ l J. aH+1 \nJl+l . f CA n+ l \ 


(y+ *x^xx- ^2 


r- 1 / .n+1 ,.n+l Ji+l »*. 


o, y 

E-2y - ^ (*f +1l j - +M,jX] + (+ 1 + 1 , J + + 1 -l,J ) [ y + ik^+i+i.j'+i-i.j]' y 


V ^2 (+u+i ‘ 2 +i!j + +i!j-i> ■ for y + +x 5 


Case 2. Shock region: ' 

[4 & * xi - < - <><] - ^x’n+r 1 - ^i 1 )- 


( y+ +x*+xx H 


0 < Xj-x 


tM E 5*1-l.j + E 6 + E 5 ( y+ ¥?-l.j +E 3+r 1 ) +E lX]^?!] 

+ (> +E 2*M,j +E 3*r 1,(E 5+(-l,j +E 6+s + 4 t 


Ax < x^- 


tn+U^n+l 


E°,( D 6+?-l ,J«7<2,j +D 8*r 1 ) +D 5' y+D 2+?-l .jXLjX' > +0 l D 5+1>U 
+ (y+D 2 +1-j,j +D 3 +^-2,J +D 4+s +1 )( [) 6+1-l,J +D 7+?-2,J +D 8+s 


2AX < X 


( 1 / ,n+l 9 ,n+l . n+1 * 

j ( +i,j+r 2 +i,j +i,j-i ) • 


"yy 


Regular point 


Ay 


. n+1 


2 ,n+1 


.n+1 


d[d+AyT C ‘ W *V . J + WfifT *U-i y J * y s < y J+i * d ' y s " y l 

liy(iWT+U+l ' 4y <} + ¥llw) C’ • y J-i < y s < y j • e “ y J " y ‘ 


Here a * x< - x s , the D n , E n are coefficients determined by appropriate Taylor series expansions. 

The A are obtained from the appropriate choice of the above representations, 
n 


+ +x J 0 
+ A a 0 

y X 

+ +x < 0 


s ‘ « 
x $ < 2 ax 
;-X s * 3 AX 


. s* 



Table 2. Illustration of Shock Fitting 
Position of shock point b: 


x b - x a 
y b ' y a 


-{ 


y + 2 (4> x + 4> x ) 


Properties ahead of the shock: 

r 

% (y+$ x ) 1/2 + (y+<t> x ) 1/2 

x be x bd 

6 = 4 - (y+d> V' 2 (<J> -6 ) 

\ % x bd X b X d 

Properties behind the shock: 

(y + 4’„) 1/,2 (^ ) = <t\, " 


) 1/2 <L + (y+4>„) 1/2 <f>„ + <t\ 


r v / \ t Y • v 

x bf x b f 


($ y 

K \ V 




— (A +d> ) L = (d) -<f) ) 

2 x b % i y b V 


V" \ '2 ( V x h )1 \\ 1 


d g 

f 


\ 

% 

\\ 

h 

A \ 


/ \ 
/ 1 

e/_ 

; 


Shock point c: 


♦c = *b + r d * * *b + tcJVV + +y b (y c- y b ) . 


Here ( ) ab means j [( ) a + ( ) b L 


Table 3. Computation Time and Rate of Convergence 


i x j 

Scheme 

Number of * 
Computations 

irk 

Time/Comp. 

Maximum 

Error 


First order 

60 

4.11 (sec. ) 

0.039593+ 

51 x 71 

Second order 

60 

4.51 

0,060134+ 


Present 

20 

9.16 

0.010475 


A computation is one complete calculation of the "solution 11 . 

■•k 

Compiling time excluded. 

Maximum error still fluctuates after 60 computations. 
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